/**** Decompose.h - Basic declarations ****/
#ifndef _H_Decompose
#define _H_Decompose
typedef struct {float x, y, z, w;} Quat; /* Quaternion */
enum QuatPart {X, Y, Z, W};
typedef Quat HVect; /* Homogeneous 3D vector */
typedef float HMatrix[4][4]; /* Right-handed, for column vectors */
typedef struct {
    HVect t;	/* Translation components */
    Quat  q;	/* Essential rotation	  */
    Quat  u;	/* Stretch rotation	  */
    HVect k;	/* Stretch factors	  */
    float f;	/* Sign of determinant	  */
} AffineParts;
float polar_decomp(HMatrix M, HMatrix Q, HMatrix S);
HVect spect_decomp(HMatrix S, HMatrix U);
Quat snuggle(Quat q, HVect *k);
void decomp_affine(HMatrix A, AffineParts *parts);
void invert_affine(AffineParts *parts, AffineParts *inverse);
void mat_mult(HMatrix A, HMatrix B, HMatrix AB);
float vdot(float *va, float *vb);
void vcross(float *va, float *vb, float *v);
void adjoint_transpose(HMatrix M, HMatrix MadjT);
Quat Qt_(float x, float y, float z, float w);
Quat Qt_Conj(Quat q);
Quat Qt_Mul(Quat qL, Quat qR);
Quat Qt_Scale(Quat q, float w);
Quat Qt_FromMatrix(HMatrix mat);
void Qt_FromMatrixNoReturn(HMatrix mat, Quat* retval);
float mat_norm(HMatrix M, int tpose);
float norm_inf(HMatrix M);
float norm_one(HMatrix M);
int find_max_col(HMatrix M);
void make_reflector(float *v, float *u);
void reflect_cols(HMatrix M, float *u);
void reflect_rows(HMatrix M, float *u);
void do_rank1(HMatrix M, HMatrix Q);
void do_rank2(HMatrix M, HMatrix MadjT, HMatrix Q);


#endif
